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We provide a scheme for efficient simulation of a broad class of quantum optics experiments. 
Our efficient simulation extends the continuous variable Gottesman-Knill theorem to a lage class of 
non-Gaussian mixed states, thereby identifying that these non-Gaussian states are not an enabling 
resource for exponential quantum speed-up. Our resuls also provide an operationally motivated 
interpretation of negativity as non-classicality. We apply our scheme to the case of noisy single- 
photon-added-thermal-states to show that this class admits states with positive Wigner function 
but negative P-function that are not useful resource states for quantum computation. 



I. INTRODUCTION 

There have been a variety of approaches to the problem of characterizing what is non-classical about quantum 
theory. There are many important signatures of quantum theory, but with the rise of quantum information, the 
exponential speedups in quantum algorithms over the best known classical algorithms have increasingly become an 
important signature of quantum theory. This point is especially poignant in light of recent work by Aaronson and 
Arkhipov wherein a simple non-universal linear optical system is shown to be able to perform computational tasks 
believed to be hard for classical computers pQ. The extent to which computational speedups and the boundaries 
between computational complexity classes are reflected in more traditional measures of non-classicality remains, 
however, an open question. In continuous variable quantum theory, and quantum optics in particular, the most 
frequently considered notions of quantumness are phrased in terms of the so-called quasi-probability distributions, 
such as the Wigner function and the (Glauber-Sudarshan) P-function. There is a strong tradition in physics of 
considering negativity of the quasi-probability function as an indicator of non-classicality of a quantum state [2H5]- 
It is therefore natural to suspect that negativity is intimately linked to computational speedups in both discrete and 
continuous quantum information processing. 

Continuous variable quantum information theory provides a potentially powerful alternative to the usual discrete 
formalism and many of the seminal results in discrete variable quantum computation have analogs in the continuous 
variable setting. Perhaps the most important example is the "continuous variable Gottesman-Knill theorem" , which 
states that a computation restricted to the subset of quantum theory containing only Gaussian states and operations 
is classically efficiently simulatable [BJ [7J . More concretely, unitary Gaussian quantum information is defined to be 
the following set of operators (see, for example, [S]): n mode Gaussian input state; quadratic Hamiltonians; and, 
measurements with (or without) post-selection onto Gaussian states. Bartlett et al. 6, 7 have shown explicitly that 
there exists a classical algorithm that reproduces the output probabilities of the measurement results that executes in 
time that scales polynomially with the number of modes. This shows that some non-Gaussian resources are necessary 
to obtain an exponential speed-up with quantum optical experiments, but leaves open the question of whether they 
are sufficient. 

Recently, Veitch et al. [S] have shown that a discrete analog of the Wigner function [TU] can be used to define 
a necessary condition for a mixed quantum state to enable an exponential speed-up through quantum computation. 
Their model considers the use of Clifford operations on audits (the case of qubits is not covered by their proof) and 
measurements of stabilizer states and finds, somewhat surprisingly, that there exist a class of bound universal states 
outside of the convex set of stabilizer states that can still be efficiently simulated and therefore do not serve as a 
resource for exponential speed-up with quantum computation. There is a tight mathematical correspondence between 
the discrete and continuous Wigner representations, the Clifford/stabilizer model for qudits 11 and the Gaussian 
model for quantum optics considered by Bartlett et al. [BJ [7] . It is therefore natural to ask whether the restriction to 
Gaussian states in the model of Bartlett et al. can be relaxed to allow more general class of initial states that have 
non-negative Wigner representation while still permitting an efficient classical simulation. 

This work affirms an answer in the positive by showing that a large class of quantum states with positive Wigner 
representation exists outside the convex hull of the n-mode Gaussian states that can be efficiently simulated using a 
classical computer, given restrictions to quadratic Hamiltonians and Gaussian measurement. This shows, in particular, 
that linear optical quantum devices are essentially no more computationally powerful than classical computers under 
such restrictions. We show this by providing an explicit classical simulation algorithm that can be used to simulate 
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sampling the output probability distributions of the evolved initial states. As a practical application we apply our 
results to determine a threshold on the computational power of single-photon-added-thermal states (SPATS) [TST 
[Ti] for variable efficiencies. In this sense, our work serves as both a conceptual and practical generalization of the 
continuous variable Gottesman-Knill theorem to a broader class of input states. 

This paper is outlined as follows. We begin with a brief review of the Wigner function formalism and Gaussian 
quantum mechanics in section |H} W e then provide our simulation protocol for states with positive Wigner represen- 
tation in section [m} In section |IV] we discuss positivity of the Wigner function as a necessary condition for quantum 
computation. We illustrate the bound state region via the recently studied class of limited-efficiency SPATS and show 
that quantum efficiencies of 50% are a necessary threshold for computational speed-up. Finally, section [V] contains 
our conclusion and further discussion about our findings. 



II. REVIEW OF WIGNER FUNCTIONS 



Wigner functions provide a natural quantum analog of the classical phase space distribution of a dynamical system. 
We provide below a brief review of the properties of Wigner functions. For simplicity, we focus our attention on 
Wigner functions for a single particle (or equivalently a single mode) in one dimension. The generalization to higher 
dimensions and more particles is straightforward [15] . 

The Wigner representation of a state p is defined to be [T5"] 

/•oo 

W p (q,p) = (27T)- 1 / (q-y/2\p\q + y/2)e^y/ h dy, (1) 



where \q) is a position eigenstate. The Wigner function is both positive and negative in general. However, it otherwise 
has many of the same properties as a classical probability density on phase space. For these reasons, the Wigner 
function is often referred to as a quasi-probability function. Intuitively, if we could find a bona fide joint probability 
distribution of non-commuting observables, then there would be no difference between quantum and classical theories. 
It is not surprising, then, that negativity is necessary in all possible quasi-probability representations of a quantum 
state [T6] . 

The time-evolution of the Wigner function for a Hamiltonian of the form H = p 2 /2m + V(q) is given by [THl HZj 

dW p (g,p) _r TTWX ,^r 1 (J\ U 9 2l+1 V(q) d 2e+1 W p (g,p) 

Ot ~ { ' pi + f^(2£+l)\\ 2) dqU+i cV^ 1 ' () 

where {•, •} is the Poisson bracket, which governs classical Hamiltonian equations of motions. This result is important 
because it states that the time-evolution of W p {q,p) is given by Liouville's equation, plus a quantum correction. The 
quantum correction is zero for the case of the quadratic Hamiltonian: 

™&>>={H,W P }. (3) 

Hence the evolution equation agrees precisely with the classical predictions. This observation will be vital for our 
simulation algorithm because the Hamiltonians permitted by linear optics are quadratic (Harmonic oscillators) , which 
(along with our non-negativity assumption) implies that we will be able to simulate the evolution of W p (q,p) using 
an ensemble of classical trajectories. 

This discussion above implies that a Wigner function that is initially classical, meaning that it is non-negative 
and hence interpretable as a probability density function (pdf), will remain classical under the action of a quadratic 
Hamiltonian. In this context is therefore useful to determine the conditions under which a Wigner function is non- 
negative as this gives a practically relevant boundary between quantum and classical states. Hudson's theorem [T%] 
was the first attempt to characterize the positive Wigner functions and it was later generalized to the following [19j . 
Let tp be a pure quantum state of n oscillators (modes) . Then its Wigner function is positive if and only if 

V-(Q) = e -h(Q-AQ+B-Q+c)^ (4) 

where A is an n x n Hermitian matrix, B is an n-dimensional complex vector and c is a normalization constant. In 
quantum optics terminology, these are either coherent states or squeezed states. That is, plugging these states into the 
definition of the Wigner function yields multivariate Gaussian distributions in phase space. Convex combinations of 
these states (incoherent mixtures of them) also have positive Wigner function since the mapping is linear. Early on, 
these were incorrectly conjectured to be the only such mixed states with positive Wigner function. The question of 



3 



mixed states was given a full treatment in reference [5D] and later in [5T]. Both references independently found that 
a theorem in classical probability attributed to Bochner [22] and generalization thereof can be used to characterize 
both the valid Wigner functions and the subset of positive ones. What was shown is that there exist a large class of 
states with positive Wigner function that are not convex combinations of Gaussian states. So far, these states have 



received little attention. In Section IV we show that such states are more than a mathematical curiosity; they arise 
naturally in quantum optics. 

Gaussian measurements are also easily modeled in the Wigner representation. Recall that for non-negative states 
the Wigner function picture allows us to represent the system as a probability density over underlying physical states 
in phase space, Uf. Gaussian measurements in this picture are also modeled as probability densities for outcomes fc, 
conditioned on the value of the underlying physical state u f . Specifically, consider the case of measurement M. of a 
Gaussian state G with covariance matrix Vm- The POVM representation of this measurement is |23) 

M(V M ) = {\G(k,V M )(G(k,VM)\ ■ fceK 2 "}, (5) 

which selects a Gaussian state with mean k and covariance Vm from all possible Gaussian states with mean k and 
covariance Vj^\ . In the Wigner function picture the representation of this measurement is, 

M M(V M )( k \ u f) = — exp (— (fe - u f ) T V^-{k - u f )) , (6) 

where X is the normalization constant. We introduce the notation in ^ to emphasize the difference between the repre- 
sentations of measurements and states. The interpretation of this equation is that if the system is actually at the point 
Uf the effect of a measurement will be to produce an outcome k according to the probability density M^/ VM ^(k\uf). 
Using this equation and the law of total probability we can find the probability density of measurement outcomes k 
for the measurement M{Vm) on a quantum state with Wigner representation W p (k): 

P (k\M{V M ),p)= [ W p (u)M M[VM) (k\u)du. (7) 

J U 

Of course, this agrees with the probability assigned by the Born rule 

p(k\M(V M ) lP ) = Tr(\G(k,V M )(G(k,V M )\ P ). (8) 

The simulation algorithm that we propose uses none of the special properties of Gaussian measurements other than 
the fact that they have a positive Wigner representation and that we can efficiently draw samples from a Gaussian 
distribution. This means that our results will apply to any measurements that satisfy these properties. We focus our 
attention on Gaussian measurements rather than these more general measurements because of their simplicity and 
physical relevance. 



III. SIMULATION ALGORITHM 



At first glance it may seem that a simulation algorithm for linear optics may be difficult owing solely to exponential 
size of the Hilbert space dimension that is generated by the evolution. We overcome this problem by exploiting the 
fact that our Hamiltonians are quadratic in p and q, which implies that the evolution of the Wigner function follows 
the Liouville equation as shown in Eq. ([3]). Since the Liouville equation preserves non-negativity and probability mass, 
the Wigner function will remain a classical distribution throughout the evolution. This allows us to model evolution 
of the Wigner function using an ensemble of classical trajectories, each of which can be efficiently simulated. The 
resulting trajectories can be used to efficiently draw samples from the final distribution of measurement outcomes 
prescribed by the Born rule without needing to know the final quantum state. 

It is critical to understand that we are not simulating the evolution of the quantum state, rather we are simulating 
measurement outcomes from a quantum circuit; this is exactly analogous to the difference between knowing a prob- 
ability distribution and being able to sample from it. In particular, the ability to efficiently draw samples does not 
imply the ability to efficiently learn the underlying distribution because the dimension of the probability distribution 
on n modes is exponentially large. 

We show in this section that this simulation strategy can be used to efficiently sample from the output of the 
following class of quantum algorithms: 
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Algorithm Class 1 Family of efficiently simulatable quantum algorithms 

Input: Number of modes n, an initial n mode quantum state p — p\ ® ■ ■ ■ <g> p n where each pj has positive Wigner 
representation W Pj (u), a description of a linear optical transformation Ut, x which is parameterized by T g fl| 2 «x2n 
and x e R 2n . 

Output: A string of measurement outcomes k sampled according to the probability density p(-K" qua nt ~ k) determined 
by the Born rule. 

1. Apply the linear optical transformation Ut,oi- 

2. Perform the separable Gaussian measurement A4(Vm) = -M(Vm,i) ® A4(Vm,2) ® ■ ■ ■ ® M(Vm,u), where we follow the 
naming convention of equation [5] and the tensor product is understood to mean that the POVM elements of M(Vm) = 
•M(Vm,i) ® • • • ® M(Vm,u) are tensor product combinations of the POVM elements of M(Vm) in the obvious way. 

3. Return the measurement outcome fc = (fci, fc2, • • • , fc„) € R 2n corresponding to the mean of a Gaussian POVM element. 



Here we conceive of any quantum algorithm in this class as a way of sampling outcome strings k distributed according 
to the probability densities given by the Born rule. We label the corresponding random variable it quant . Here we are 
not simulating the evolution of the Wigner distribution, which would be equivalent to simulating the full quantum 
state. Rather, we show that there is a corresponding classical algorithm that produces outcome strings k with (very 
nearly) the same distribution those from algorithm class [I] 

Using intuition similar to that in [5] , we note that quantum algorithms in class [l] can be simulated using the 
following classical algorithm, provided access to classical resources with infinite numerical precision and a blackbox 
function that can be used to draw samples from W Pj (it) for j = 1, . . . , n. We will later provide an algorithm that does 
not require infinite precision, but we provide the infinite precision algorithm first because it conveys the necessary 
intuition without focusing on the technical issues that arise when discretizing the distributions. 

Algorithm 1 Infinite precision classical simulation algorithm for algorithms in class [T] 
Input: As algorithms in class 111 except p is not provided. 

— (2) 

Output: A string of measurement outcomes k sampled according to the probability density p(K^ ss = k). 

1. Sample u £ R 2 ™ according to the distribution W p (u) — W pi (ui) ■ ■ ■ W Pn (u n ) by drawing a sample independently from 
each mode using the blackbox function. 

2. Apply the afhne transformation u = Tu + x corresponding to the linear optical transformation Ut,hs to the sampled 
phase space point u. This transformation is an affine mapping due to Louiville's theorem. 

3. Return the outcome string fc = (fci, fc2, • ■ • , k n ) G K 2n from the distribution 

M M(v M )(k\u) = M m{Vm il )(fci|i*i) ■ • ■ M M ( VM n) (k n \u n ), where M M (v M>i )(ki\ui) is given as in equation^ 

M Mv Mij ){ki\ui) = exp (-(ki - u t ) T Vy^j(fci - it;)) . 



The intuition behind this class of algorithms is to use the classical phase space model afforded to us by the 
non-negative Wigner functions and quadratic evolutions in order to turn the quantum problem into one that can be 
efficiently simulated by a classical computer. In this context we can think of our quantum system as actually being 
definitely at some point it £ M. 2n which is unknown to us. The point then moves under a fully deterministic classical 
evolution and measurement on each register amounts to picking a point k from a normal distribution centered at the 
location of the system. Each classical algorithm samples outcomes k according the probability density 

P(K%L s = k)= I W p (u)M M{VM) (k\u)du, 

J u 

which agrees with the density given by the Born rule. Thus the outcomes -Kpiass fr° m the classical simulator are 
distributed in exactly the same way as -Kquant, which are outcomes drawn from the actual quantum system. 

Unfortunately, algorithms similar to [l] cannot be executed precisely on digital computers and instead would require 
an analog computer (often referred to as a real computer). If physical, such computers would have unrealistic 
computational powers such as being able to solve NP-complete problems in polynomial time [24] and would also 
violate the holographic principle [25] , For these reasons, we need to discretize[l]in order to assess the cost of simulating 
linear optics on realistic classical computers. The major technical difference between the continuous variable case and 
the discrete case considered in [!|| involves showing that finite-precision errors can be made negligible with efficient 
overhead costs given a set of reasonable assumptions about the input states, dynamics and measurements. 
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Since infinite precision is requried in the continuous variable setting, in order to specify a quantum state we must 
make some finite precision truncation. To this end we shall assume access to a family of oracles W Pi >ri (l, m) that takes 
integers l,m and returns a value satisfying \W Pjt ^(l,m) — W Pi {n Pj + (l,m)S)\ < r\. That is, each oracle queries the 
Wigner function at points on a grid centered at the mean of the distribution. This is a weak assumption as it does 
not require us to even know the state we are simulating. Using this resource the algorithm can be written as: 



Algorithm 2 Finite-precision classical simulation algorithm for quantum algorithms in class [T] 

Input: As algorithm [l] but also require S, a discretization length for the input, | €2 1 , a bound for the numerical error 
involved in applying the affine transformation, T, a discretization length for the output, \A\, the area truncated square 
region of phase space that the simulator considers and \i Pj , the mean of the Wigner distribution of the quantum state 

on each mode j. We require \J\A\ to be an odd integer multiple of 5 and Y to be an odd integer multiple of 5. 
Output: A string of measurement outcomes k sampled according to Pr(i^ class = k) = Pr sim (fc). 

1. For each j = 1, . . . ,n execute a through d. (This step approximates sampling a point from phase space.) 

set Pr sim , Pj (I, m) = W Pj , v (l, m) ■ S 2 . 

(The phase space is truncated to a size |.4| and discretized into boxes of size 8. This step sets a pdf over the centers 
of the boxes.) 

(b) For each (l,m) set Pr S i m , p . (I, m) = Pr s i m ,p . (I, m)/ Y^i m P r sim,p,- {I, tti) (This step normalizes the pdf.) 

(c) Draw a sample (l,m) from the pdf Pr s i mjPj . (I, m). 

(d) Set Uj = fij + (I, rn)5 

2. Set u = Tu + x using enough digits of precision such that the numerical error is at most | e2 1 , where u = ui © • • • u n 
and u = Ui © • • • © U n - 

(This step corresponds to updating the sampled state according to the linear optical transformation.) 

3. For each j — 1, . . . ,n execute a through e. (This step is to simulate drawing a measurement outcome from the Gaussian 
measurement distribution centered at u.) 

set Pi sim ,M(v M ,j)( l ' m ) : = ex P {- s (h , m) T Vj^j5(l,m)))) ■ 5 2 . 

(The outcome space is truncated to a size |^4| and discretized into boxes of size 8. This sets a pdf over the centers 
of the boxes.) 

(b) Set PT slmiMi v M<i ){l,m) = PT aimtM( y M)j )(l,m)/Y^i <m ^uim,M(v M>j )ihm). 
(This step normalizes the pdf.) 

(c) Draw a sample (l,m) from the pdf Pr s i m . .m(v_m )(l> m )- 

(d) Find integers (r,s) such that \5(l, m) - T(r, s))|oo < T/2. 
(This just amounts to rebinning the outcome distribution into hypercubes of sidelength T; this introduces no errors 
but does require 8 < T.) 

(e) Set measurement outcome kj = iij + T(r, s). 

4. Return measurement outcome fe = ki © • • • © fe n . 



Our simulation protocol can necessarily only sample from a discrete distribution so we must introduce some notion 
of how a discrete distribution can be close to the continuous probability density given by the Born rule. The most 
natural way to do this is to discretize the outcome distribution into boxes of side length T according to, 

Definition 1. Let Afr(k) C M 2 " be a hypercube in outcome space with side length T centered at the point k. 
We define the T discretization of the quantum outcome distribution to be Pr(-Kq U ant,r = k) = Pr qua nt,Ai(yM)(' i; ) = 

Xvr(fc) P(-^q uant = fydk. 

We can now fix a Y according to our operational requirements for the simulation and ask how well a simulation pro- 
tocol samples from this distribution. Notice this is an unavoidable consequence of trying to approximate a continuous 
quantity with a discrete system. With this in hand we can give a precise classical simulation protocol by discretizing 
our naive algorithm, which results in the family of protocols described in algorithm [2j 

It now easy to see that both the cost of the simulation and the error in our sampling will be a function of the 
discretization parameters S(n, e) and |-4(n, e)|. If we can pick these parameters such that for fixed error our simulation 
scheme scales as poly(n) then the simulation is efficient. Concretely, 



(a) For each integer I. 



m £ 



(a) For each integer l,m £ 



G 



Definition 2. Let Pi S iin.,M(v M ) (&) an d P r quant,A-i(y.M) C 5 ) ^ e simulated and actual probabilities of obtaining a 
measured value that is inside a hypercube of volume T 2n centered at a point in phase space k £ K 2n for the separable 
Gaussian measurement M(Vm)- A simulation algorithm is efficient if for inputs n, T and e there exists a choice of 
\A\, 6 such that: 

1. The 1-norm distance between the T-discretized quantum distribution and the simulator distribution is at most 

I " r quant , M ( Vm ) 

(k)\i<e 

where we take PT S i m ,M(v M )(k) = whenever k is outside the domain of PTsim,M(v M )(k) defined by algorithm[2j 

2. The computational complexity of the simulation scales as poly(n/e). 

We now can show that the simulation of sufficiently smooth separable positive Wigner functions under linear optical 
operations and Gaussian measurements is efficient. This result is formally stated in the following theorem and proof 
is given in appendix [A"| 

Theorem 1. The output of algorithm [2] satisfies |Pr quant m(v m ) (&) — P r sim M(v M )(tyh — e f° r input 
({V Pj },{V Pl },{VMA,T,x,6,\e 2 \,P,\A\Yii 

1. n, ma,Xj{\(i Pj |, \\V Pj , \\VM,j\\}, \\T\\ and \\x\\ are bounded, 

2. There exist finite /3, A such that \VW p (u)\ < n/3/\A\ n a,nd \W k M M(VM) (k\u)\ < nA/\A\ n for k £ R 2n , uel 2 ", 

3 - ^ min { i 6 [(i + ||T||)A + ^^ ' r }' 

4. e<l, |e 2 | < \\T\\8y/*, 

5. \A\ > lGnmax.j ([V Pi } 11 + [V Pi ] 22 + [V M ,j] u + [V M , 1 } 22 ) / e , 

6. The finite precision error from each oracle yV Pj . v satisfies rj < gra |^i n . 

Furthermore, if we assign unit cost to evaluations of W Pj tV and Gaussian functions unit cost then the computational 
complexity of the algorithm is O (n 5 (max, ||VpJ] + max^ 1 1 V^vi , j 1 1 ) [A 2 ||T|| 2 + /3 2 ] /e 3 ) , which implies efficiency. 

The key insight of this theorem is that the assumption of Gaussian preparations made in the continuous variable 
Gottesman-Knill theorem can be relaxed [51 [7] , and that a much wider class of quantum dynamics can be efficiently 
simulated than previously thought. Indeed, although we stated the algorithm only for product state inputs and 
product measurements we can now see that this restriction was unnecessary. In fact, our simulation scheme works for 
any positively represented input and measurement as long as it is possible to efficiently sample from the corresponding 
distribution. The product assumption is a sufficient but not necessary condition for this efficient sampling. We also 
note that the algorithm requires us to know the mean and covariance matrices of the distributions, which might 
be hard to compute analytically. However, since we already require efficient sampling we may appeal to Monte 
Carlo estimation protocols to determine these quantities within acceptable error tolerances. This extension of the 
continuous variable Gottesman-Knill theorem places much stronger limitations on the input states that can be used 
for continuous variable quantum computation and underscores the significance of negativity in the Wigner function 
as a resource for quantum computation (in analogy to recent results for discrete systems). In particular, we will show 
that Theorem [l] places a minimum efficiency for a class of photonic thermal states beyond which the states cannot be 
used as a resource for linear optical quantum computation with Gaussian measurements. 



IV. EFFICIENT SIMULATION OF SINGLE PHOTON-ADDED-THERMAL STATES 

The debate over the "correct" definition of classicality for quantum states of light has been a long and, at times, fierce 
one. One of the most common notions of classicality is whether a state can be represented as a convex combination 
of Gaussian states. Here we have shown that, in the context of computational power, such a condition is superseded 
by the condition of positivity of the Wigner function. In this section we give a concrete example of an interesting 
class of states which are not Gaussian but which have positive Wigner representation and thereby admit an efficient 
classical simulation. 
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We consider the experimentally accessible class of states called single-photon-added thermal states (SPATS) [T/THT1] : 

^ oo 



nmm, 



where n is the mean photon number — given by, in terms of temperature T, the Planck distribution n = 1/ (exp(l/T) — 
1). It is known that all states in this class are outside the convex hull of Gaussian states and have negative Wigner 
function for finite temperatures. 

Under experimentally realistic conditions we must consider states subjected to losses. In general, losses can be 
modeled as an interaction with a vacuum mode at a beam-splitter with transmittance rj, also called the quantum 
efficiency. The loss rate is then 1 — r\. The Wigner function of the SPATS after this channel, which we call limited- 
efficiency SPATS (or LESPATS for short) is [26] 

, , 2 1 + 2r][n + 2(n + \){q 2 + p 2 ) - 2 w? - 1] / 2(q 2 + p 2 )\ 
W »^P) = ~ TT^-p 6XP tlWj- 

Note that the most negative value of the LESPATS is at (q,p) = (0, 0) for all r\ and n. Thus, we consider the quantity 

2 1 + 2r?(n - 2nrj - 1) 



tt (1 + 2m?) 3 



By inspection, we can see that for efficiencies of ij < 0.5, the Wigner function of the LESPATS is positive W p ( HiT/ ) (<?, p) > 
0. Thus, efficiencies of rj > 0.5 are necessary for quantum computational speed-up with SPATS. 

Note however that for rj < 0.5, the LESPATS are not inside the convex hull of Gaussian states. To see this most 
clearly, we require a different phase space distribution. The Glauber P-function (see for example [27 ) is defined as 

P p {q,p)\a)(a\dqdp, 

where \a) are the coherent states, which are vacuum states that have been displaced in phase space (symmetric 
Gaussian states). Note that if P p is a probability distribution then p is in the convex hull of coherent states. The 
P-function of the LESPATS is [H] 



irn r) 



(n+1) 



exp 



q 2 +p 2 
fvq 



Note that this function is negative for all allowed values of ij and n. Thus, the LESPATS are always outsides the 
convex hull of Gaussian states but are bound universal states [9] for r\ < 0.5. 

To illustrate this, we compare the negativity of the Wigner function with the distance to the convex hull of Gaussian 
states. The distance we use is based on fidelity, which can be computed using phase space distributions as 



F(p(n,r,), |0)(0|) = Tr[p(n, r,)\0)(0\] = tt J J P p{%v )(<l,p)Q\o)(o\(q,p)dqdp, 

where the Q-function 

Q\o)(o\(q,P) = ^cxp(-(<? 2 +p 2 )) 

is dual to the P-functiorQ Using a standard table of Gaussian integrals we find 

1 — T) 



F{p{n,rf), 



(l + my) 2 



This effect is demonstrated in figure [T] Note that, for any state, a quantum efficiency of 77 < 0.5 is sufficient to 
ensure membership of the convex hull of states that have positive Wigner representation. This effect is mirrored in 
the discrete case [9] , where depolarizing noise of 50% is sufficient to ensure membership of the polytope of states with 
positive discrete Wigner function, when starting from any qudit state. 



1 This duality relationship holds for any phase space representation 1281 
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Figure 1. Negativity of Wigner function on the left and fidelity to vacuum ("distance" to convex hull of Gaussians) on the 
right for varying n and n (note: n — 1 corresponds to no losses). In both figures "+" indicates the region of non-negative 
states and "-" indicates the region of states with negative Wigner function. The region of non-negative states (n < 0.5) is the 
region of bound universal states. This is clear as the Wigner function is positive yet the states still lie outside the convex hull 
of coherent states since the P-function is always negative. Notice that the fidelity distance to the convex hull (as measured by 
the fidelity to the nearest state, |0)) is significantly less than 1, suggesting that the region of bound states is quite large. 



V. CONCLUSION 



We have shown that Gaussian quantum computations utilizing separable initial preparations with positive Wigner 
function are classically efficiently simulable. Since such states lie outside the convex hull of Gaussian states, we have 
identified a large class of bound states: states that cannot be prepared using Gaussian operations, yet do not permit 
universal quantum computation. We illustrated this class using the example of single-photon-added-thermal-states, 
showing that quantum efficiencies of 50% are necessary for quantum computation. 

Effort has been extended beyond qualitatively defining negativity as quantumness to quantifying quantumness via 
negativity. In terms of the Wigner function, the volume of the negative parts of the represented quantum state 
has been suggested as the appropriate measure of quantumness [29] . The distance (in some some preferred norm 
on the space of Hermitian operators) to the convex subset of positive Wigner functions was suggested to quantify 
quantumness in reference [30]. The volume of negativity of the Wigner function (and, in the finite dimensional case, 
the sum of the negative values) is a "good" measure of non-classicality since it is monotonic under Gaussian operators; 
that is, Gaussian operations cannot increase the volume of the negative regions in phase space. 

Reference [31] nicely summarized what was known at the time about continuous variable quantum computation. 
The table presented there is reproduced below in table [I] with some more recent results. The field began with 
Lloyd and Braunstein's observation that non-linear optical processes are sufficient for universal continuous variable 
quantum computation. Later, it was shown for discrete variable encodings that linear optics is sufficient provided 
photon counting measurements are available [32, 33 . The continuous variable analog of the measurement-based model 
shows that preparation of single photon state preparation is also sufficient [M] . More recently, the result of Aaronson 
and Arkhipov [1 shows that preparing and measuring single photon states (without post-selection) is equivalent to 
a sampling problem that is thought to be hard classically — but it still manages to (probably) not be universal for 
quantum computation. 

It is possible that the Aaronson and Arkhipov model may be intermediately between classically efficiently simulat- 
able and universal for quantum computation. Another suspected model of this type is the "one-clean-qubit" model 
of Knill and Laflamme 35 . The key point for this latter model is that uses highly mixed states. Mixed states have 
not been given full consideration for continuous variable quantum computation. Here we have shown, via the Wigner 
phase space formalism and independent of purity, negative representation is necessary for universal quantum com- 
putation. Moreover, any computation that uses states possessing a positive Wigner function is classically efficiently 
simulatable. It would be quite interesting if this condition turned out also to be sufficient as this would provide a 
sharp boundary between quantum and classical systems with regard to their computational power. 
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Appendix A: Proof of Theorem \T\ 

Proof of Theorem^ Our goal is to show that with this choice of discretization parameters S and \A\ the simulation 
algorithms outlined in algorithm [2] require 0(poly(n/e)) resources and result in error at most e. Since we bin the 
data at the end of the protocol into hypercubes of volume T 2n > 5 2n and T is an integer multiple of S, no error is 
introduced by first binning the outcomes into hypercubes of volume S 2n because every hypercube of volume S 2n is 
contained in precisely one hypercube of volume T 2n . We therefore may take the quantum distribution to be binned 
into hypercubes of side length 5 (ie. r = S) without loss of generality. 

We start by analyzing the error. Following scheme outlined above we can decompose this into four broad parts: 

1. The error introduced by the use of finite precision output of W Pji?) . 

2. The error introduced by truncating the sampling distribution over phase space. 

3. The error introduced by discretizing this truncated distribution. 

4. The error introduced by truncating the outcome distribution. 

5. The error introduced by discretizing the outcome distribution. 

Denoting the region in phase space that the initial states are confined to as A p — A Pl ® • • • ® A Pn and the region 
that the observations are confined to as Am{v m ) = •A-m(v m ,-l) ® " ' ® •^■M(y M , n ) ( wnere \A Pi \ — \AM{v M ,i)\ = 
we can use the triangle inequality to express these errors as: 



\Pr S im,M(V M ){k) — Pr q uant,7Vl(y M )(fc)|l < | Pr S imJV( (V M ) (&) ~ Pl'trunc-quaiit,A1(y. M ) (&) 1 1 

+ Prquant^ £ A p ) + Pr quant (fc £ A M (V M )), (Al) 

which says that the total error is at most the distance between the truncated distributions plus the probability 
that a point is sampled, or measured, outside the truncated region. In other words, the 1-norm error introduced 
by truncating is, even in the most pathological case conceivable, the sum of the probabilities of sampling an initial 
trajectory outside of A p and measuring an outcome outside of Am(v m A- 
It then follows that \Pr silnM ( VM) (k) - Pr quant ,M(VCM)( fe )li ^ e / 2 if 

Pr qU a„t(M i A p ) + Pr quant (fc £ A m[Vm) ) < e/4, (A2) 

\^sim,M{V M )( k ) ~ P r trunc-quant,A^(V M )(fc)| 1 < e/4. (A3) 



We will first demonstrate that ( A2 ) is an immediate consequence of assumption [5] We then will show that ( |A3[ ) is 
satisfied if assumption [3] holds. 

We begin by bounding the truncation error. For the i th mode, 



/ \/\A\ \/\A\ 

Prquant («i £ A pi ) = Pr |% - fJ, Q \ > , \ Pl - fj,p\ > 



< Pr \ qi - MQ | > ^ + Pr \ Pi - fi P \ > 





Our upper bounds for both of these probabilities are established using Chebyshev's inequality Chebyshev's inequality 
states that, for a probability distribution P with mean \x and standard deviation a that P(\x — fx\ > ka) < (a/k) 2 . 
In our case, the two standard deviations are -\/[Vp;]ii and -\/[Vpi] 22 which gives us 

Pr qU a nt (^^ P J< 4([ ^ ]l ; + [ ^ ]22) . (A4) 
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It immediately follows from the independence of the n distributions that compose W p that, 

[V Pi ] n + [V Pi } 22 



By identical reasoning, 



Pr q uant(M A p ) < in max 

i— l,...,n 



Piquant (fc ^ A M (V M )) < 4n max 

j = l,...,n 



Thus we have, 



Pr qU ant(fc i A p ) + Pr quant (tt A M {V M )) - 4nmax 



1-41 

[^M,j]ll + [Vm,j}22 

\A\ 

[V Pi ] n + [V Pi ) 22 + [V M ,j]ii + [Vm.M 



\A\ 



so choosing | ^4. | > 8nm&Xi 



[^] 11 + [ y P,] 22 + ^^hl + ["KM, 3 ]22 



e/2 



guarantees, 



Pr q uant(M 4- Ap) + Pr qU ant(fc ^ A M (V M )) < 



(A5) 



We must now bound the discretization error on the truncated distributions. It is natural to break this error up into 
three pieces corresponding to the first three steps of the algorithm. First, there is the error introduced by discretizing 
the initial sampling distribution over phase space. Next there is the numerical error introduced in implementing the 
affinc transformation u = Tu + x. Finally, there is the error introduced by discretizing the outcome distribution. The 
remainder of the proof is devoted to bounding these three errors and thereby bounding the total error using ( Al I. 



Step 1 of the simulation algorithm breaks up A p into hypercubes of volume 5 and samples from the set of centers 
of these hypercubes according to, 



Pr s im, P (M = M + (h, 



, l n , mi, . . . , m n )S) = Pr siniiPl (l u mi) • • • Pr sim , p „ (l n , m„) 
= Wp 1>17 (Zl,7Tli) ■ • • Wp niV (l n , m n ) ■ 5 2n 
= W pl (u 1 )---W Pn (u n )-S 2n + e 1 , 



(A6) 



where Ci is the numerical error introduced by using n > 0. The probability weighting assigned to a hypercube center 
u is only approximately equivalent to the probability mass contained in the hypercube; we must bound the error 
introduced by this nonequi valence. Define Ns(u) C K 2n to be the hypercube with side length 8 centered at the point 
u. The for a fixed u in the domain of Pr s i m p (tt) we have that every point w € J\f$(u) is also in the region of truncation 
A hence W p :r j unc (w) — W Pj (w) for all such w. We then use this simplifying observation, the mean value theorem and 
the triangle inequality to find, 



Pr. 



sim,p 



(«)- / 



W p (v)dv 



< 



< 



max \VW p (v)\ 



max 

uGAf s (u 



\v — u\ 



r2n+l 




61. 



(A7) 

\VW p (v)\) < n/3/\A\ n by assump- 



Where max^gTv^^jIt; — u\ < \/2n^ from Pythagoras' theorem and (max^g^^) 
tion [2] of Theorem [Tl 

In the second step of the simulation algorithm we move the sampled point u £ M. 2n to the point u = Tu + x + £2 , 
simulating the evolution due to the linear optical transformation Ut, x - The numerical error e 2 depends only on 
numerical precision which can be made exponentially small using a linear amount of memory. The other source of 
error in the simulation of the dynamics is due to error in the initial conditions caused by sampling the point u S M. 2n 
as opposed to the point v £ M. 2n that would have been sampled if i5 — 0. The triangle inequality then implies that 



(T« + a; + |e 2 |)-(r^ + a;)| < |e 2 



T|| max \v 

vEATs(u) 



£2 



imia 



(A8) 



Since e 2 can be made exponentially small using a polynomial number of computational steps, we can choose | e2 1 < 
||T||<5a/¥ without affecting the efficiency of the algorithm. Therefore, by making such a choice, the total error in the 



simulated dynamics is at most 



51 < \\T\\5V2n~. 



(A9) 
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The final step of the simulation algorithm samples a measurement outcome on each mode. To do so we break up 
the truncated outcome space Am(v m ) into hypercubes of volume S 2n and sample from the set of centers of these 
hypercubes according to, 



Pr s im,A^(Vx)( fe = A + ( l U ■ ■ ■ ,ln,mi, . 



,m n )6\u) = Pr sim ,x(y Mil )(ii,mi) • ■•Pr slm ,M(VM, n )(^ 

_ 1 
_ l 



i exp (- J2 5 ( l i, m l ) T V J ^ i 5{k, m^j • S 2 
^ exp (- jhfa - ty T Vj$y M j(ki - Ui)J ■ 5 2 



(A10) 



where C is a normalizing constant and fci are components of k. As in step one we may use the mean value theorem 
to bound the error introduced by sampling this way rather than according to the true probability mass over each 
hypercube. Using Afs(k) C K 2 ™ to be the hypercube with side length 5 centered at the point k we find from the 
mean value theorem and the assumptions of theorem [l] that the error in the probability enclosed in a single hypercube 
centered at k is at most: 



Pr 



sim,VW (Km) 



(k\u) 



M M(V M ){lAu)dK 



< max 



|V ' K M M (y M ) 



(k\u)S 



2n I 



max \k - k\ 

KS7V 5 (fe) 



< 



h5 2n+1 ny/2n 



\A\» 



(All) 



We complete the error analysis by bounding the distance between the simulator distribution and truncated, dis- 
cretized quantum distribution. Our aim is to rewrite this error in terms of quantities we have already bounded. Let 
T> u be the domain of Pt^ m j^(v M )(u); i e - the discrete set of points in A p that can be sampled by the simulator. 
The difference between the probabilities of a fixed outcome k occurring for both the simulation and the truncated 
quantum distribution 



|P r sim,M(VM)(fc) — P r trunc-quantJV!(yM)(&)| 

Y] PTshn,M(V M )(k\u + e 2 )Pr simtP (u) - / / 

JAfs(k) JA P 



M^ c Mi) (K\v)W p (v)dvdK 



(fc|u + e2 )Pr sim ,»- Y, I I Mj^ Mi) (K\v)W p (v)dvdK 



< Y Pr sim,M(V M )(k\u 



£2) 



(u)- / W p (v)dv 



Y] / Pr s im,M(v M )(k\u + e 2 ) - / 



W p {v)dv. 



(A12) 



The final inequality is found by adding and subtracting X^ex>„ P^sim.M(v M ) (k\u + e 2 ) /^(u) W p {v)dv and applying 
the triangle inequality. 

We need one more intermediary bound before arriving at the final result. Let v £ J\f$(u). It follows from the mean 
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value theorem and M^ u ,?? 



L 



M s (k) 



M_m(v m ) throughout the domain ol integration that 
Mm&m) + e ^ dK - I M^ M) {K\v)dK 



\ M m(Vm) (k\u + e 2 )dK - / M M (V M ){n\v)dK 



< 



< 



JVs(fe) 



< 6 



max 
nA 



M s {k) \A n 
2n nA 



V u M^ m) (k\u) 



\T(u -v) + e 2 \dn 



Ms{k) 

\T(u ~v)+e 2 \dn 



\T\\S 



\T\\S 



(A13) 



where we have used the assumption that \e 2 \ < ||T||(5-y/n/2 and identical reasoning to bound \T(u — v)\. In conjunction 
with (All I this implies: 



P^sim, M(v M ){k\u + £2) - / 

J A 



<(i+imi)^^ — +n 



(A14) 



Since there are numerical errors in the computation of the probability distribution for measurement outcomes, it 
is conceivable that the algorithm assigns more than unit probability to outcomes in the region A (although this is 
unlikely for small values of e). Although we can claim that J^uev Xv 5 (u) W p (v)dv < 1 we can only claim that 
Suex> ^ >I sim,M(v M )(k\'U' + e 2) 5: 1 + e < 2 under the assumption that e < 1. Using these facts, we substitute (A7) 
and AAlty into (|A12| to find: 



\P*sim,M(V M )(k) -Prtrunc-quant,A<(V M )(fe)| < K 1 + 1 1 T \ I ) A + I 3 } 



nV2nS 2n+1 



\A\» 



Since we have assumed that the outcomes are discretized into hypercubes of volume 5 2n , the discretization produced 
^L- hypercubes. This implies that the 1-norm distance between the the two distributions is: 



rune— quant, J\4 {Vm ) (fe)| < [(l + ||T||)A + /3]n\/2^+| ei ||.4|7<5 2n 

k 

Since there are n modes and numerical error n, it is straightforward to see that |ei| < nr]S 2n ; therefore choosing 

c 



S < 



and 



implies 



8[(1 + \\T\\)A + f3]nV2n' 



rj < 



8n\A\ r 



|P r sim,.M(V.M) (fc) P r trunc— quant, M{V M ) (fc)ll < 



(A15) 



This sets the discretization error to be e/4. We found previously that the truncation error is at most e/4 under 
assumption [5] of theorem [I] We then combine ( |Al[ ), (A5) and ( A15| ) and find: 

|P-^sim. M(v M ){k) - Pr 

quant, M{Vm) 

Cfc)U <e/2, 

if r = S. This also implies that the result holds for T > 5 given that T is an integer multiple of 5 (which ensures that 
the hypercube of size S 2n that the measurement is assigned to is inside the hypercube of size T 2n that it should be 
assigned to given T-discretization); therefore, it generally holds if 



S < 



8[(l + ||T||)A + /3] nV2n 



(Ai6) 
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A final issue remains: we did not consider errors that are introduced in the sampling steps in the algorithm that 
arise because the simulated probability distributions are not normalized. We will show that, from the assumption that 
e < 1, it suffices to divide e by 2 in all previous calculations. To see this, let y be a discrete probability distribution 
and let s be a vector such that \y — x\% < e' for some 1 > e' > 0. It then follows from the triangle inequality that 
1 + e' > \xU >l-e'. Thus 



\y - x /\ x \i\ 1 = \y\ x \i - x \ 



= \y\x\i — x\x\i + x\x\i — x 
< \y-x\i + \\xU -II < 2e', 



(A17) 



under the assumption that e' < 1. Therefore, by combining these results with (A16) we see that the difference between 
the (now normalized) simulated probabilities and the quantum predictions is at most e if 

S < min | — t= , T \ , (A18) 



16[(l + ||T||)A + /3] W2^' 



and 



\A\ > 16nmax ([V Pi ] n + [V p .} 22 + [V M j] u + [V M ,j} 22 ) A, (A19) 

as claimed by theorem [l] 

We complete the proof by showing that algorithm [2] is computationally efficient with this choice of \A\ and S. Again 
we will analyze the simulation algorithms step by step. In the first step we sample a point u t on the phase space 
of each register from the distribution Pr S i m p (ttj). This distribution has support on ^ squares, and thus if we take 
\A\ and S to be proportional to their respective lower and upper bounds then the number of times that W PjtV must 
be evaluated scales as (n 4 (||Vp|| + ||Vm||) [A 2 ||T|| 2 + (3 2 ] /e 3 ) (here 0(-) is Bachmann-Landau notation meaning 
asymptotically bounded above and below by a constant multiplied by (•)). If we ascribe unit computational cost to 
every such access, then the computational complexity of this step is proportional to the number of times that W PJlV 
is queried. This task must be repeated n times, and hence the total computational complexity of this step is. 

e(n 5 (||^|| + ||yM||)[A 2 ||T|| 2 + /3 2 ] a 3 ). 

In step two we apply the affine transformation u — ¥ Tu + x, whose cost is dominated by the cost of performing a 
matrix multiplication using 0(log(||T||i^/n)) bits of precision. Since the matrix multiplication requires a number of 
arithmetic operations that scales quadratically with the matrix dimension and addition and multiplication scale at 
most quadratically with the number of bits of precision, the total cost of this step is 0(n 2 log 2 (||T|j<5y / n)), which is 
subdominant to the cost of the previous step and therefore does not affect the scaling. 

In step three, we are confronted with the task of measuring the resultant trajectory using a separable Gaussian 
measurement. There is a strong duality between drawing a sample trajectory from the initial separable Wigner 
function and drawing a measurement outcome for the separable Gaussian measurement of the final trajectory. In 
particular, we discretize the space surrounding the outcome space of each of the n separable measurements into |-4|/<5 2 
points. The approximate calculation of the measurement probability requires that we perform a number of operations 
that are proportional to the number of points. Therefore, identically to step one, the computational cost is 

0(^ 5 (II^II + IIKm||)[A 2 ||T|| 2 + /3 2 ] /e 3 ), 

which verifies the computational complexity claimed by theorem [l] and shows that under the assumptions of the 
theorem all three steps are computationally efficient, and hence the algorithm is efficient as well. 
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Preparations 


Gates 


Measurement 


Efficiently simulatable 
classically 


Vacua 


Linear optics 


Gaussian 


/ [6] [7] 


Vacua 


Non-linear optics 


Gaussian 


X [37] 


Single photons 


Linear optics (no 
squeezing) 


Photon counting (with 
post-selection) 




Vacua 


Linear optics 


Gaussian and Photon 
counting (with post- 
selection) 


X [33] 


Single photons 


Linear optics 


Gaussian 


X [38] 


Single photons 


Linear optics (no 
squeezing) 


Photon counting 




Product Positive Wigner 
functions 


Linear optics 


Product Gaussian 


/ (this work) 











Table I. An extension of the table appearing in [31] . 



